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Abstract We simulate the rise through the upper convection zone and émer- 
gence through the solar surface of initially uniform, untwisted, horizontal mag- 
netic flux with the same entropy as the non-magnetic plasma that is advected 
into a domain 48 Mm wide from from 20 Mm deep. The magnetic field is advected 
upward by the diverging upflows and pulled down in the downdrafts, which 
produces a hierarchy of loop like structures of increasingly smaller scale as the 
surface is approached. There are significant différences between the behavior of 
fields of 10 kG and 20 or 40 kG strength at 20 Mm depth. The 10 kG fields have 
little effect on the convective fiows and show little magnetic buoyancy efîects, 
reaching the surface in the typical fluid rise time from 20 Mm depth of 32 hours. 
20 and 40 kG fields significantly modify the convective flows, leading to long 
thin cells of ascending fluid aligned with the magnetic field and their magnetic 
buoyancy makes them rise to the surface faster than the fluid rise time. The 20 
kG field produces a large scale magnetic loop that as it émerges through the 
surface leads to the formation of a bipolar pore-like structure. 

Keywords: Sun; Dynamics; Magnetic fields; Magnetic Flux Emergence; Con- 
vection 



1. Introduction 

We have initiated simulations of the rise of magnetic flux through the near 
surface layers of the solar convection zone. Numerous studies exist of the rise 
of a thin magnetic flux tube through the bulk of the convection zone, e.g. 
dMoreno-Insertis, 1986; Choudhuri and Gilman, 1987"; Fan, Fisher, and McClymont, 1994fl 
[Moreno-Insertis, Caligari, and Schuessler, 1994; Caligari, Moreno-Insertis, and Schussler, 1995] ) .[ 
Thèse breakdown near the bottom of our simulation because the magnetic pres- 
sure exceeds the gas pressure and the tubes explode. Others have studied the 
rise of cohérent, twisted, non-thin flux tubes through model convection zones, 
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e.g. dDorch, 2QQ3||Abbett et al, 2QQ4|rFan, 2QQ9D . The rise of a cohérent, twisted 



flux tube through the shahow layers of reahstic solar convection zones has been 
studied by (|Cheung et a/., 200 8"; Cheung, Schûssler, and Moreno-Inserti s, 2QQ7| 



Cheung, Moreno-Insertis, and Schûssler, 2006; Martmez-Sykora, Hansteen, and Carlsson, 2QQ9H 



Martmez-Sykora, Hansteen, and Carlsson, 2008). We are undertaking a corn- 



phmentary study - the rise and évolution of untwisted, horizontal flux with 
the same entropy as the unmagnetized ascending plasma, carried by upflows 
from a depth of 20 Mm. Here we présent some preliminary results from thèse 
simulations. 



2. Simulation 

Our domain is 48 Mm wide and from the température minimum to a depth 
of 20 Mm (Fig. [T]). We solve the équations for mass, momentum and inter- 
nai energy in conservative form plus the induction équation for the magnetic 
field, for fully compressible flow, in three dimensions, on a staggered mesh. 
The code uses finite différences, with 6th order derivative operators and 5th 
order interpolation operators. The grid is uniform in horizontal directions with 
resolution 95.24 km and non- uniform in the vertical (stratified) direction with 
resolution of 12 km near the surface increasing to 75 km at the bottom of the 
domain. Transformation between uniform index space, where the derivatives 
are evaluated, and the non-uniform physical grid is achieved by multiplying by 
the Jacobian of the transformation. The staggered mesh increases the order of 
the derivative operators for the same size stencil, but the physical variables are 
not aligned in space and need to be interpolated using fifth-order interpolation 
operators. Time intégration is by a 3rd order low memory Runge-Kutta scheme 
( [Kennedy, Carpenter, and Lewis, 1999D . Parallelization is achieved with MPI, 
communicating the three overlap zones that are needed in the 6th and 5th order 
derivative and interpolation stencils. Because of the staggered mesh, div(B) is 
conserved to machine précision. However, numerical errors in div(B) that slowly 
accumulate stochastically must be cleaned at intervais. Ghost zones are loaded at 
the top and bottom boundaries to permit the use of the same spatial derivative 
scheme at the boundaries as in the interior. Horizontal boundary conditions are 
periodic, while top and bottom boundary conditions transmit the convective 
flows but partially reflect waves. Inflows at the bottom boundary have uniform 
pressure, specified entropy and damped horizontal velocities. Outflow boundary 
values are obtained by extrapolation. The magnetic field is made to tend toward a 
potential field at the top. That is, the ghost zones are loaded with a potential field 
calculated from the values of the vertical field just inside the boundary. At the 
bottom the magnetic field is given a specified value in inflows and extrapolated 
in outflows. The electric field in the ghost zones is calculated from the magnetic 
field and the velocities (which are taken to be constant at their values on the 
boundary at the top and have zéro vertical derivative at the bottom) in the 
ghost zones. The magnetic field is updated from its time derivative, which is the 
curl of the electric field. The code is stabilized by diffusion in the momentum, 
energy and induction équations. 
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Figure 1. Mean atmosphère of the model. The domain includes ail of the hydrogen and heli- 
uml ionization zones and most of the hélium II ionization zone (H dashed, Hel dash-dot, Hell 
dash-dot-dot-dot, Fi solid). It covers only 10% of the géométrie depth of the solar convection 
zone, but half of the pressure scale heights. 




Ionization energy accounts for 2/3 of the energy transportée! near the solar 
surface and must be included to obtain the observed solar velocities and tempéra- 
ture fluctuations ( |Stein and Nordlund, 1998| ). We use a tabular équation of state, 
that includes local thermodynamic equilibrium (LTE) ionization of the abundant 
éléments as well as hydrogen molécule formation, to obtain the pressure and 
température as a function of log density and internai energy per unit mass. 

Radiation from the solar surface produces the low entropy, high density fluid 
whose buoyancy work drives the convective motions and it also provides us 
with our only information on what is occurring there. Since the optical depth 
is of order unity in thèse régions, neither the diffusion nor optically thin ap- 
proximations are valid. We calculate the radiative heat/cooling by solving the 
radiation transfer équation in both continua and lines using the Feautrier method 
dFeautrier, 1964D , assuming Local Thermodynamic Equilibrium (LTE). The num-| 
ber of wavelengths for which the transfer équation is solved is drastically reduced 
by using a multi-group method whereby the opacity at each wavelength is placed 
into one of four bins according to its magnitude and the source function is binned 
the same way dNordlund, Ï982l|Stein and Nordlund, 2QQ3|rVQgler, Bruis, and Schûssler, 2QQ4D .| 
The conversion of energy transport from convection to radiation occurs in a 
ver y thin thermal boundary layer, which must be resolved (grid size less than 
15 km) in order to obtain the correct entropy jump and the solar convective 
adiabat. This radiation solution is validated by comparing a few snapshots with 
a detailed solution having a large number of wavelength points and many more 
angular quadrature points. Calculated spectral line widths, Doppler shifts, and 
asymmetries from the simulations are in excellent agreement with observations 
dAsplund et al, 2QQQD . 

At super granule and larger scales, the coriolis force of the solar rotation begins 
to influence the plasma motions, so /-plane rotation is included in the simulation. 
Angular momentum conservation produces a surface shear layer with the surface 
(top of the domain) rotating slower than the bottom of the domain, as observed 
in the Sun. 

We do not know the properties of the magnetic fleld at the depth of 20 Mm. 
We have investigated the behavior of untwisted, horizontal flux with the same 
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Figure 2. Vertical fîeld (bottom row) and émergent continuum intensity (top row) at times 
33, 40, and 48 hours after the 10 kG field started to be introduced at the bottom. The range 
in the field strength is [-2.8,2.84] kG. 

entropy as the unmagnetized ascending plasma, advected into the computational 
domain by upflows from a depth of 20 Mm. Wherever there was inflow at the 
bottom, uniform horizontal field in the x-direction is advected into the computa- 
tional domain. We have calculated three cases, with field strength at the bottom 
of 10, 20 and 40 kG. In the outflows the field values are extrapolated. At the 
boundaries between in- and out- flows, the field magnitude is smoothly varied. 
The boundary conditions are applied to the electric field from which the time 
derivative of the magnetic field is calculated. The initial state in ail three cases 
was the same snapshot of field free hydrodynamic convection. 

3. Flux Emergence 

Upflows carry the field toward the surface, and for the case of the stronger fields, 
magnetic buoyancy increases the speed of ascent. The typical time for fluid to 
rise from a depth of 20 Mm to the surface is 32 hours. The field in the stronger, 20 
kG, case reaches the surface in about 20 hours, assisted by magnetic buoyancy. 
The field in the weaker, 10 kG, case reaches the surface in the fluid travel time 
of 32 hours. Simultaneously, convective downflows pull portions of the field lines 
downward. Turbulence twists and stretches the fieldlines. As the magnetic flux 
approaches the surface it spreads out. In the weaker (10 kG) case, the downflows 
are able to confine the field to only a portion of the 48 Mm wide domain, but 
in the stronger (20 and 40 kG) cases the field spreads to over the entire domain 
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Figure 3. Vertical field (bottom row) and émergent continuum intensity (top row) at times 
19, 20 and 23 hours after the 20 kG field started to be introduced at the bottom. The range 
in the field strength is [-3,3.4] kG. 



at the surface. It is quickly swept by the diverging upflows to the boundaries of 
the granules and more slowly to the boundaries of meso-granule scale structures 
(Fig. [2]). The emerging magnetic field distorts the granules, elongating them in 
the direction of the horizontal field. As the large loop in the 20 kG case reaches 
the surface, it produces a large bipolar région (Fig.[3|) with a typical field strength 
of 2 kG at continuum optical depth unity. The simulations have not yet run long 
enough to see if the supergranule scale diverging flows will produce a magnetic 
network. 

At the bottom of the domain, the horizontal magnetic field is entering the do- 
main in the supergranule scale upflows. Thèse upflows, plus magnetic buoyancy, 
cause the field to rise to the surface. The downflows at thèse early times keep 
the field pinned down near the bottom of the domain. As a resuit, the downflows 
are nearly field free and the magnetic field, even though untwisted, is confined to 
the upflows (Fig.[4j). As the surface is approached, the vertical velocity becomes 
broken up into smaller and smaller cells and the magnetic field is spread around. 
However, the area above the large deep downflow still remains mostly field free. 
It becomes ringed with field lines diverging from it on one side and converging 
toward it on the other (Fig. [5]). This is the manifestation of the very large loop 
structure that is produced (Fig. [8|). While the weaker field is rising, the large 
scale downflow pattern at depth is changing and the magnetic field does not 
have such cohérent large scale loops. 
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Figure 4. Image of vertical velocity (red and yellow downward, blue and green upward) and 
horizontal magnetic field during the initial rise through the non-magnetic atmosphère for the 
20 kG case at depth 8 Mm and 4 Mm. At thèse intermediate depths the magnetic field is being 
advected in the upfiows. The downfiows are still confining the magnetic fiux to deeper layers 
and so are mostly field free at thèse depths. The field is slowly spreading horizontally with 
decreasing depth. 




Figure 5. Image of vertical velocity (red and yellow downward, blue and green upward) and 
horizontal magnetic field during the initial rise through the non-magnetic atmosphère for the 
20 kG case at depth 2 Mm and 1 Mm. Nearer the surface the magnetic field has spread farther 
and now almost surrounds the area of the large downfiow at deeper levels. 



4. Magnetic Topology 

The magnetic field is transported toward the surface in the upfiows (assisted by 
magnetic buoyancy where the field strength is large enough) and pushed down 
in the downdrafts which créâtes Q and U loops from the initially horizontal field. 
Since the size of the horizontal cellular structure of convection flows decreases 
as the surface is approached (from supergranule size near the bottom to granule 
size near the surface, Fig. [6]) this produces loop structures over this range of 
scales as well (Fig. [7]). 
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Figure 6. Vertical velocity's horizontal pattern at the surface and depths 2, 4, 8, 12, and 16 
Mm below the surface. Blue and green are upflows. Red and yellow are downflows. The rms 
velocities are respectively 1.47, 0.93, 0.55, 0.35, 0.31 and 0.24 km/s. The size of the horizontal 
cellular patterns of the vertical velocity increases with increasing depth as the scale height 
increases. 




Figure 7. Example magnetic field line for the 10 kG case, showing loop like structure on 
many scales, from 2 to 20 Mm. Volume visualization is scaled vertical velocity (to make its 
visibility uniform with depth). Blue and green upwards, yellow and red downwards. Distance 
scales are in Mm. 
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Figure 8. Large scale magnetic loop extending the full domain in the x-direction (direction of 
the horizontal field entering the domain from below) and from the bottom at 20 Mm to almost 
the surface. The volume rendering is the scaled vertical velocity (to make the flow at depth as 
visible as at the surface). Red and yellow are downflows and green and blue are upflows. The 
axis labels are megameters. 




Figure 9. Small magnetic loops near the surface. The plane near the top is the visible surface 
with the vertical velocity pattern. Note that the field lines (yellow, green, red and blue) that 
are looping over a granule interior connect very differently beneath the surface. 



An example of a large scale loop structure spanning most of the 48 Mm width 
and 20 Mm depth is shown in Fig. [8l The magnetic field was buoyant as well 
as being transported by the upflows and was pushed down in the downflows. 
An example of small loops near the surface is shown in Fig. [9l The field lines 
close to the surface are very twisted on the scale of a few megameters. Note that 
the several field lines that loop over a granule interior connect very differently 
beneath the surface because of the turbulent convective motions. 

In the 20 kG case, there is enough magnetic flux that it becomes concentrated 
in large régions to field strengths of 2 kG at continuum optical depth unity and 
produces pores (Fig.[Tn]). The intensity in the pores is 1/4 of the average intensity. 
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Figure 10. Emergent continuum intensity (left) and vertical magnetic fîeld at continuum 
optical depth unity (right) in a plage région of the 20 kG case. 




Figure 11. Magnetic Field lines and émergent intensity in bipolar pore région looked at from 
above and below. In this case many field lines are connected in a U-configuration. 



The magnetic field lines in this pore région are generally nearly vertical at the 
surface and have a U-shape below the surface (Fig. [TT]). The opposite polarity 
field lines at the surface generally connect at large depths (greater than 10 Mm) 
below the surface. 

In ail three cases the magnetic energy density is larger than the fluid kinetic 
energy density below 10 Mm. However, this is not necessarily an argument 
against such strong fields, since they could be confined because they are in pres- 
sure equilibrium with their surroundings. Thèse preliminary resuit s, however, 
suggest that the magnetic field can not be as strong as 20 kG at 20 Mm unless it 
has a lower entropy then the field free plasma, because magnetic bubbles develop, 
due to the large magnetic buoyancy, and produce large, hot, bright granules at 
the surface, which are not seen. More cases need to be calculated and compared 
with observations to clarify the properties of the magnetic field at 20 Mm depth. 

5. Application to Helioseismology 

/— and p— mode waves are excited in the simulation and can be used for testing 
and refining local helioseismic techniques. Figure [T2l shows an example of waves 
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Figure 12. Acoustic waves (visualized by the density fluctuation they produce) in the non- 
magnetic simulation propagating into the interior and refracted by the increasing sound speed 
with depth as calculated by J. Zhao. 




0.0 2.0 4.0 6.0 8.0 0.0 2.0 4.0 6.0 8.0 

k (Mm"') k (Mm"') 



Figure 13. k — uo diagram comparing non-magnetic and magnetic (20 kG) case. On the left 
the magnetic contours are overlaid on the non-magnetic image and vice- versa on the right. 
There is no shift in the mode ridges and the image on the left looks like the contours on the 
right and vice versa, so that there does not appear to be any différences produced so far in 
this calculation. 

gêner ated near the surface, propagating into the inerior, and refracted by the 
increasing sound speed with depth. We have made prehminary k — uo diagrams 
and compared the magnetic and non-magnetic cases. For the early stage of the 
flux émergence studied so far there are no discernable différences fFig. [T3|). 

6. Discussion 

The results we have reported here concern only the intial rise of uniform, un- 
twisted, horizontal magnetic flux through the near surface, non-magnetic solar 
convection zone. The main resuit, so far, is that fields as strong as 20 kG at 20 
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Mm depth are not possible unless they have lower entropy than the surrounding 
field free plasma. Using thèse simulations, we will continue to investigate how 
the magnetic flux evolves over time, how the pores evolve, whether sunspots 
form, whether a magnetic network is produced and what its properties are, how 
flux émerges through an existing magneto-convecting surface layer, how thèse 
simulations compare with observations and what constraints can be placed on 
the properties of the magnetic field near depths of 20 Mm. The /— and p— 
modes excited in thèse simulations will be used to investigate local helioseismic 
inversion procédures. 
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